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Abstract Solar energetic particles (SEPs) are one of the extreme space weather 
phenomena. A huge SEP event increases the radiation dose received by aircrews, 
who should be warned of such events as early as possible. We developed a warn¬ 
ing system for aviation exposure to SEPs. This article describes one component 
of the system, which calculates the temporal evolution of the SEP intensity and 
the spectrum immediately outside the terrestrial magnetosphere. To achieve this, 
we performed numerical simulations of SEP transport in interplanetary space, in 
which interplanetary SEP transport is described by the focused transport equa¬ 
tion. We developed a new simulation code to solve the equation using a set of 
stochastic differential equations. In the code, the focused transport equation is 
expressed in a magnetic field line coordinate system, which is a non-orthogonal 
curvilinear coordinate system. An inverse Gaussian distribution is employed as the 
injection profile of SEPs at an inner boundary located near the Sun. We applied 
the simulation to observed SEP events as a validation test. The results show that 
our simulation can closely reproduce observational data for the temporal evolution 
of particle intensity. By employing the code, we developed the WArning System 
for Aviation Exposure to Solar energetic particles (WASAVIES). 
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1 Introduction 

Space radiation poses a serious threat to several human activities, such as high 
radiation doses for astronauts, adverse effects on aircrew health, artificial satellite 
malfunctions, and the disruption of high-frequency radio wave communications in 
high-latitude regions. Predicting the space radiation environment to reduce the 
risk of radiation hazards is one of the most important goals of space weather 
research. 

Space radiation is primarily composed of energetic protons and electrons. Ra¬ 
diation doses from the protons are harmful for astronauts and aircrews, whereas 
the electrons predominantly affect artihcial satellites. Energetic protons originate 
from galactic cosmic rays (GCRs) and solar energetic particles (SEPs). The spec¬ 
trum and intensity of GCRs, and hence their radiation doses, are almost constant 
over time scales of much less than the 11 year solar cycle. On the other hand, 
SEPs occur sporadically and their intensity suddenly increases by multiple orders 
of magnitude within days. Such sudden increases in SEP intensity present a major 
radiation hazard to astronauts, adversely affecting the success of space missions 
[Hu et al., 2009]. While low-energy SEPs cannot penetrate a spacecraft or ex¬ 
travehicular activity suits, high-energy SEPs, the energies of which exceed several 
tens of MeV, penetrate these materials [Kronenberg and Cucinotta, 2012; Reames, 
2013]. Therefore, predicting high-energy SEP events is the first priority in ensuring 
the safety of astronauts and the success of space missions. For an excellent review 
on the space radiation environment, refer to Vainio et al. J2009]. 

Despite the lack of understanding of the SEP mechanism, many researchers 
have attempted to predict the occurrence of SEP. However, a physics-based dehni- 
tive prediction remains difficult because SEP occurrence is related to various phys¬ 
ical processes yet to be fully clarihed, such as solar flares, coronal mass ejections, 
coronal and interplanetary shock waves, and particle acceleration and transport 
mechanisms. Therefore, previous studies on predicting the occurrence of SEP have 
adopted statistical, empirical, or probabilistic approaches Je.g., Garcia, 2004a, 
2004b; Kubo and Akioka, 2004; Kuwabara et al., 2006; Kahler et al., 2007; Posner, 
2007; Balch, 2008; Laurenza et al., 2009; Nunez, 2011]. 

On the other hand, Reames J2004] stated that, “In the author’s opinion, reliable 
predictions of the onset and fluence of an SEP event prior to its occurrence are not 
likely in our lifetime. However, after an event onset, it should be possible to model 
the CME, the shock, and the acceleration and transport of particles sufficiently 
well to predict the peak intensity at shock passage and the event fluence.” Some 
researchers have focused on simulating the temporal evolution of SEP intensity, 
and numerous simulation studies of particle acceleration and transport in inter¬ 
planetary space have been published Je.g., Lario et al., 1998; Zank et al, 2000; 
Rice et al., 2003; Sokolov et al., 2004; Verkhoglyadova et al., 2009]. However, the 
aim of most of these studies was not to predict the SEP intensity prohle. To the 
best of our knowledge, the only operational prediction model that incorporates the 
physical mechanism of SEPs was developed by Aran et al. [2006], who predicted 
the SEP intensity profile from numerical simulations of interplanetary shock prop¬ 
agation and energetic particle transport. Their model is specific to low-energy SEP 
intensity profiles, the energy of which is less than several tens of MeV. 

A ground level enhancement (GLE), which is one of the extreme space weather 
phenomena, is induced by extremely energetic SEPs having energies greater than 
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450 MeV (approximately 1 GV in rigidity for protons) [Shea and Smart, 2012]. 
Energetic SEPs can produce secondary particles such as neutrons, muons, neu¬ 
trinos, and gamma rays through nuclear reactions in the terrestrial atmosphere. 
All such secondary particles induce a marked increase in the radiation dose rate 
at aviation altitudes. This means that effort should be made to alert aircrews to 
the presence of extremely energetic SEPs so that they can take necessary action 
to reduce the radiation dose rate they are exposed to. Recently, Kataoka et al. 
[2014] developed the WArning System for AVIation Exposure to Solar energetic 
particles (WASAVIES). The system aims to warn aircrews of high dose levels and 
to simulate GLEs by forward modeling from SEP transport to an air shower in 
the terrestrial atmosphere [Sato et al, 2014]. 

This article describes a newly developed numerical simulation code to calculate 
SEP transport in interplanetary space for the WASAVIES system. The next section 
describes the SEP transport equation and a coordinate system to express the 
equation, along with its solution method. Then, the results of test calculations are 
described. Subsequently, comparisons of calculations with observational data are 
presented, which are followed by a discussion section. The final section presents a 
summary. 


2 Numerical simulation of solar energetic particle transport 

SEPs are believed to be accelerated within three regions: the reconnection region 
in solar flares, coronal shock waves, and interplanetary shock waves. According 
to Caprioli and Spitkovsky [2014], the exponential cutoff energy of an ion energy 
spectrum produced by diffusive shock acceleration, which roughly represents the 
maximum energy achievable, is proportional to the background magnetic field in¬ 
tensity. Thus, the majority of energetic SEPs cannot be accelerated by interplan¬ 
etary shock waves, for which the background magnetic field is weak. This means 
that they tend to be accelerated near the Sun (by solar flares and/or coronal 
shocks) and then transported to the Earth along interplanetary magnetic fields. 
Therefore, energetic SEP events will be reproduced by solving the SEP transport 
equation from the Sun to the Earth if an accelerated energetic SEP time profile 
near the Sun is assumed. In this section, we describe the numerical simulation of 
SEP transport through interplanetary space. 


2.1 Focused transport equation 

As the scattering mean free path of interplanetary SEP transport is roughly com¬ 
parable to the distance between the Sun and the Earth, SEPs are not isotropically 
distributed in interplanetary space, particularly in the initial phase of SEP events. 
Therefore, interplanetary SEP transport cannot be modeled using Parker’s trans¬ 
port equation, which assumes an isotropic particle distribution. Instead, the SEP 
transport equation must explicitly incorporate changes in pitch-angle. For this rea¬ 
son, a focused transport equation is widely used to simulate interplanetary SEP 
transport [e.g., Ruffolo, 1995; Zhang et al., 2009; Droge et al., 2010; He et al., 
2011; Qin et al., 2013]. In this study, we employ the focused transport equation, 
which is a form of the Fokker-Planck equation, described as 
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dtf + nvhdif + V^dif + pdpf + (idpf - dp {Dppdpf) = 0, (1) 

where p, v, and p are the momentum intensity, speed, and pitch-angle cosine of 
a particle, respectively, i denotes the three components of spatial coordinates, bi 
and Vi are the unit vector along the magnetic field line and solar wind velocity, 
respectively, Dpp denotes the pitch-angle scattering coefficient, d with subscripts 
t, i, p, and p denote partial differential operators with respect to time, spatial co¬ 
ordinates, momentum intensity, and pitch-angle, respectiveljQ, and p and p denote 
the time derivatives of p and p, respectively. The gyrotropic phase space density 
/ is a function of i, p, p, and t. On the left-hand side of equation (HD. the second 
to sixth terms describe particle streaming, convection induced by solar wind, adi¬ 
abatic deceleration, the change in the pitch-angle due to adiabatic focusing and 
diverging solar wind, and pitch-angle scattering, respectively. The adiabatic decel¬ 
eration and change in the pitch-angle in equation (HD are, respectively, dehned as 
[Isenberg, 1997] 

b^bJd^VJ - j , (2) 

and 

1 — ^ 

p = —[-vbidi (InB) - p {SbibjdiVj - d^V^)] , (3) 

where B is the magnetic held intensity. The hrst term in the square brackets on the 
right-hand side of equation (ED describes the adiabatic focusing effect. In this set 
of equations, the particle position and momentum are measured in the co-rotating 
frame of the Sun and the solar wind frame, respectively. 

We employ a simple Parker spiral magnetic held with radial solar wind speed 
Vr and angular speed of solar rotation 17 for an interplanetary magnetic held 
structure dehned as 


p = p 


1 - Sp^ 


B = Bs (sr - , (4) 

where <P denotes the angle between the radial and magnetic held directions and is 
expressed as 

tan^ = arsin0, (5) 

where a is the ratio of 17 to V7, r, 6, and p are, respectively, the radial distance, 
polar angle, and phase angle in spherical coordinates, and are unit vectors 
denoting the radial and phase angle directions, respectively, and Be is the radial 
component of the magnetic held at a specihc distance re- 

The pitch-angle scattering coefficient Dpp is dehned by a quasi-linear theory 
[Jokipii, 1966] with some modihcation to avoid zero scattering at /r = 0 [Beeck and 
Wibberenz, 1986; Bieber et al., 1994] and is described as 

Dpp = Dovp'^~'^ (^\p\'^~^ + hj (^1 - p'^^ , ( 6 ) 

^ In a mathematically precise description, the partial differential operators with respect to 
spatial coordinates should be written as covariant differential operators, and covariant and con- 
travariant vectors should be distinctively written in a curvilinear coordinate system. However, 
these differences in description are not distinguished in this article. 
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where Do is related to the intensity of the turbulent magnetic field and q is the 
inertial-range spectral index of the tnrbnlent magnetic field. Here, a Kolmogorov 
spectrnm is assnmed, i.e., q = 5/3, and the modification constant h is set to 0.2 
[Qin et al, 2006]. The pitch-angle scattering coefficient is related to the parallel 
mean free path afl 


Am = 


8 i_i 


2\2 


3n (1 - M ) 


D, 


-dji. 


( 7 ) 


The radial mean free path is defined as 


Arr = A|| cos^ (8) 

which is assumed to be constant across the entire interplanetary space [Bieher et 
al., 1994]. 


2.2 Coordinate system 

As a convenient coordinate system for this problem, a magnetic field line coordi¬ 
nate system (s, 0, </>) is defined, as shown in Figure[T] In the figure, s is the distance 
of position P measured from the origin O along a magnetic field line, 0 is the polar 
angle at P measnred from the Z axis, which is fixed at 7r/2 in this article, and 4> is 
the phase angle of the magnetic held line at O, rather than at P, measnred from 
the X axis. The relationship between the magnetic held line coordinate system (s, 
9, (/) and the spherical coordinate system (r, 9, ip) is expressed as 


r 

^=2 


1 


H—In (tan (p + 


1 


cos sin (p 


cos(p 


= 0 , 

= OLr -\- ip. 


( 9 ) 

( 10 ) 

( 11 ) 


The magnetic held line coordinate system is identical to the spherical coordinate 
system when the solar rotation vanishes, i.e., when a = 0. 

The solar wind velocity in the co-rotating frame of the Sun only has an s- 
component in the magnetic held line coordinate system, when is expressed as 


Hs = 


Vr 

COS<P 


( 12 ) 


As the magnetic held line coordinate system is a non-orthogonal curvilinear 
coordinate system, off-diagonal elements appear in the metric tensor, and the line 
element dl^ of the coordinate system is described as 


dl^ = ds^ -|- r^dO^ -|- (rsin 9)^ dcj)^ — — tan^ ^cos^ dsdcj), (13) 

a 


^ A small correction is required for the parallel mean free path defined by equation © in 
the case of focused transport [Danos et ai, 2013 and references therein]. However, equation 
0 is used as the relation between the pitch-angle scattering coefficient and the parallel mean 
free path in SEP transport simulations in this article. 
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where the relationship between r, s, and 9 is given by equation ©■ Equation o 
can be recast in the magnetic field line coordinate system as 

dtf + {^J-v + 14) dsf +pdpf + (A - dpDpp) d^f - Dp,pdp,dp,f = 0, (14) 

where the rate of change of the momentum and p and rate of change of the pitch- 
angle A are, respectively, expressed as 


Vs cos 

p = p - 


1 - 

2 


sin^ <P - 



(15) 


and 

A = —cos^ ^2 -h tan^ -|- pl4 [2 — tan^ j . (16) 

Equation (HI can easily be solved because of the spatial one-dimensionality along 
a magnetic field line. 


2.3 Method used to solve focused transport equation 

The focused transport equation described above is solved by solving a set of 
stochastic differential equations that is mathematically equivalent to the Fokker- 
Planck equation [Kloeden and Platen, 1999; 0ksendal, 1999; Zhang, 1999]. Stochas¬ 
tic differential equations are increasingly used to model phenomena such as SEP 
transport [e.g., Zhang et al., 2009; Droge et al, 2010; Qin et al., 2013], cosmic-ray 
transport in the heliosphere [e.g., Pei et al., 2010; Strauss et al., 2013], cosmic-ray 
transport in the Galaxy and interstellar media [e.g., Farahat et al., 2008; Effen- 
berger et al., 2012], and particle acceleration in shock waves [e.g., Zuo et al., 2011]. 

While the Fokker-Planck equation can generally be recast as a set of both 
time-forward and time-backward stochastic differential equations, we use a set of 
time-backward stochastic differential equations because a time-backward method 
is suitable for calculating the time evolution of particle intensity or the energy 
spectra of an SEP event observed near the Earth [Kopp et al., 2012]. Therefore, 
equation (US is recast as a set of time-backward stochastic differential equations 
as follows: 


ds = — {pv -(- 14) dr, (17) 

dp = —pdr, (18) 

dp = - {fi- dpDpp) dr + ^/TO^dw, (19) 


and 


dr = —dt, (20) 

where t is the time measured backward and dw in equation m denotes a one¬ 
dimensional Wiener process, the probability density function of which is described 
by the Gaussian function 
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p {dw- dr) = exp 


dw^ 
2 dr 


( 21 ) 


To numerically solve this set of stochastic differential equations, the motion 
of pseudo particles can be tracked by considering numerous stochastic processes 
using a Monte Carlo method. In accordance with equation m, random variables 
are generated for each step of numerical integration by a pseudo random number 
generator. Pseudo particles are numerically traced backward in time over a speci¬ 
fied duration. When a particle comes in contact with an inner or outer boundary 
of the calculation domain, calculations for that particular particle are terminated. 
The inner and outer boundaries are set at heliocentric radial distances of 0.05 and 
80 AU, respectively. 

The expectation value of the generated stochastic processes weighted by an 
inner boundary condition is given by 


/ = (/o {St-To,Pt-To,^J't-To,t - To)) , 


( 22 ) 


where the brackets (• • •) denote the expectation value and tq is the backward time 
at which the pseudo particle comes in contact with the inner boundary, /o denotes 
the inner boundary condition, which expresses the SEP injection profile near the 
Sun. The expectation value f{s,p,fi,t) describes the phase space density at time 
t of the focused transport equation, which is the final solution. 

The SEP injection profile at the inner boundary is given by 



(23) 


and it is normalized to unity over the integral with ranges pmin < p < oo, 0 < 
/r < 1, and 0 < t < oo. At the time of injection, the pitch-angle is uniformly 
distributed between 0 and tt/2, and the momentum spectrum is a power law with 
index — 7 . The time profile is expressed as an inverse Gaussian distribution with 
mean m and standard deviation a. Although most of the studies on interplanetary 
SEP transport [e.g., Droge, 2000; Qin et al., 2006; Zhang et al, 2009] use the so- 
called Reid-Axford profile [Reid, 1964] as the injection time profile, we employ the 
inverse Gaussian distribution as the injection time profile for the following reason. 
The inverse Gaussian distribution is related to the diffusion process as dX = 
Xdt+5dw, where X, t, and dw denote the state variable, time, and Wiener process, 
respectively, A is the average rate of change of the state variable, and denotes 
the diffusion coefficient. The first and second terms on the right-hand side of the 
equation are called the drift and diffusion terms, respectively. The time t at which 
X first attains state A/, starting from X = Xq at t = 0, has an inverse Gaussian 
distribution. In physical terms, if energetic particles are injected impulsively at a 
specific position in the lower corona and move outward at a constant drift rate 
accompanied by diffusion, the particles that escape at a specific distance in the 
upper corona have an inverse Gaussian distribution in time. Therefore, we adopt 
the inverse Gaussian distribution as the particle injection time profile. 

Because the pseudo particles integrated backward in time from the Earth to the 
Sun are pulled toward the outer boundary by pitch-angle focusing, most pseudo 
particles cannot return to the Sun. To allow particles to return to the Sun, a 
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method proposed by Qin et al. [2006] is used, in which the focused transport 
equation is rewritten to suppress pitch-angle focusing. The expectation value is 
then calculated as the weighted average of the stochastic processes in accordance 
with the rewritten equation. The weights are calculated while the stochastic pro¬ 
cesses are generated. Mathematically, the solutions obtained using this method 
are identical to those of the original transport equation. 


3 Test calculations 

By solving the focused transport equation described in the previous section us¬ 
ing the injection profile equation, we can simulate the SEP event intensity and 
anisotropy profile. The anisotropy A{t) parallel to the magnetic field is related to 
the pitch-angle distribution as 


A{t) 




(24) 


Figure [2] shows the simulated SEP intensity (top panel) and anisotropy (middle 
panel) for 142 MeV protons at distance of 1.0 AU from the Sun for radial mean 
free paths of 0.2, 0.8, 1.4, and 2.0 AU at 142 MeV, indicated as red, green, blue, 
and pink lines, respectively. The bottom panel shows the time profile of particle 
injection near the Sun with m = 2.5, a = 5.0 hours, and 7 = 5, which are the same 
for all mean free paths. We can see from the figure that the SEP intensity and 
anisotropy strongly depend on the radial mean free path, namely the pitch-angle 
scattering. As shown in the top panel, the intensity decreases almost exponentially 
for the smallest mean free path (0.2 AU) and deviates further from exponential 
decay as the mean free path increases. As shown in the middle panel, the longer 
the mean free path, the slower the anisotropy decreases. This phenomenon demon¬ 
strates that more particles arrive from the solar direction as the mean free path 
increases because hardly any energetic particles are scattered back from behind 
the Earth. 

Figure [3] shows the temporal evolutions of the pitch-angle distribution for mean 
free paths of 0.2 AU (top panel) and 1.4 AU (bottom panel), which are shown as 
red and blue lines in Figure [2] respectively. The black, red, green, and blue points 
in the figure show the pitch-angle distributions at 1, 3, 6, and 12 h, respectively, 
following particle injection near the Sun. The pitch-angle is measured from the 
direction of the magnetic field line away from the Sun. In both cases, the pitch- 
angle gradually becomes more isotropic from its initial anisotropic distribution. For 
the 1.4 AU mean free path, the pitch-angle remains highly anisotropic even after 
12 h, whereas that for the 0.2 AU mean free path is almost completely isotropic. 
This phenomenon is attributed to the fact that hardly any particles are scattered 
back from behind the Earth owing to the long mean free path or weak scattering 
in solar wind. 

Figure [I] shows calculated rigidity spectra at the Earth for a mean free path 
of 0.8 AU at 1 GV. The eight rigidity spectra shown in the figure are for times 
of 30, 40, 50, 60, 90, 180, 360, and 720 min after the start of particle injection 
near the Sun. The dashed line shows the injection spectrum at the Sun with a 
power-law index of -5. The injection parameters are the same as those for Figure 
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[2] It is evident that the rigidity spectra change with time. This complex behavior 
of the spectral shape is ascribed to the time-of-flight effect, which means that the 
particle fluxes in the low-rigidity range are low in the early phase because low- 
rigidity particles arrive near the Earth later than high-rigidity particles do. Our 
simulation can calculate the complex behavior of SEP momentum spectra, which 
is essential to calculate neutron monitor count rates during a GLE, for example. 


4 Comparison with observations 

In this section, we compare observational data with the results of our numerical 
simulation. The data we used are differential intensities obtained at 1-min intervals 
with the GOES P 6 and P7 channels, whose energies are 165 MeV and 433 MeV 
(580 MV and 1,000 MV, respectively in rigidity for protons). Because the GOES 
satellites are located inside the terrestrial magnetosphere, the observed SEP in¬ 
tensities may have been affected by the magnetosphere. However, this effect will 
be small at the high particle energies of the P 6 and P7 channel data. To reduce 
the effect of the longitudinal difference between the two GOES satellites on the 
SEP intensities, the data collected from GOES 13 and 15 are averaged. Noise is 
eliminated by calculating the running average of data points. 

The method of comparison is quite simple; the results of the numerical simula¬ 
tions are htted to the observed data using a selected radial mean free path. Data 
are fitted using a nonlinear least-squares method and the square error is minimized 
using the simplex method [Nelder and Mead, 1965]. The mean m and standard 
deviation a of the inverse Gaussian injection prohle are adjusted to match the 
simulated intensity of the observed SEP event. Because the results are insensitive 
to the momentum spectrum index 7 , it is fixed at 5. 


4.1 Events used for comparison 

For reproducing extreme events, we selected events with intensities for SEPs having 
energy greater than 100 MeV observed by the GOES satellite reaching 1 PFU. Nine 
such events have been observed since 2012 . By checking SDO/AIA 94 — A images, 
it was found that four of the nine events occurred in the eastern hemisphere or 
near the central meridian region of the Sun (longitude less than W30), and we 
discarded these four events. One of the remaining events had missing data for 
both GOES 13 and 15. Therefore, four events remained: those on 27th January 
2012, 13th March 2012, 17th May 2012, and 6 th January 2014. We compared the 
four events with the results of numerical simulations. 


4.2 Event on 27th January 2012 

Figure [ 5 ] shows a comparison of the SEP event that occurred on 27th January 2012 
with the results of our numerical simulation. The upper and lower panels represent 
the P6 and P7 channel data, respectively. The gray points indicate the observation 
data, which are plotted in the figure by thinning them out to increase their visibil¬ 
ity. Unfortunately, as the GOES satellites do not observe the anisotropy of SEPs, 
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we do not know the scattering mean free path in interplanetary space, and the 
comparisons shonld be performed over a range of possible radial mean free paths. 
Therefore, we compared the simnlation resnlts for the SEP intensities calculated 
using the proposed method (colored lines) for mean free paths ranging from 0.2 
to 2.0 AU (increment: 0.3 AU) with the GOES data. The figure shows excellent 
agreement between the observational data for both the P6 and P7 channels and 
our simulation for all mean free paths. 


4.3 Event on 13th March 2012 

Figure |6] shows a comparison between the event observed with the GOES P6 and 
P7 channels on 13th March 2012 and the calculated intensity. The gray points and 
colored lines have the same meanings as those in Figure O The upper panel shows 
that the P6 channel data are almost perfectly reproduced by our numerical simu¬ 
lation results. In the lower panel, the P7 channel data are again well reproduced, 
although the late phase data of the event deviate slightly from the simulation 
results. 


4.4 Event on 17th May 2012 

The SEP event that occurred on 17th May 2012 was the Erst ground level enhance¬ 
ment (GLE) within the 24th solar cycle, which was an extreme space weather event. 
Figure 0 shows a comparison of the observed data with the calculated intensity. 
The gray points and colored lines have the same meanings as those in Figure [S] 
Again, our simulation results reproduce the observed SEP event with high accu¬ 
racy for all mean free paths and for both P6 and P7 data, even though the decay 
rate of this event is considerably faster than that of the event on 27th January 
2012 . 


4.5 Event on 6th January 2014 

The SEP event that occurred on 6th January 2014 was an interesting event. No 
increase in X-ray flux was detected by GOES at the time of the event. However, a 
flare was clearly observed in the eastern hemisphere from the view of STEREO-A, 
suggesting that the SEP event was related to a flare that occurred on the far side 
of the Sun. Figure [8] shows the fitting result for this event. Excellent agreement 
between the observation data and our calculations can be observed in both the 
upper (P6 channel) and lower (P7 channel) panels. As another SEP event occurred 
toward the end of this event, the data during this period were not used in the fitting 
procedure. 

The comparisons shown in this section indicate that the inverse Gaussian time 
profile can be reasonable for the SEP injection profile near the Sun. 
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5 Discussion 

We used the inverse Gaussian distribution as the SEP injection profile in the 
solar corona in this study. As already mentioned, however, the Reid-Axford profile 
is often used for the injection profile. The time dependence of the exponential 
functions included in both the inverse Gaussian and Reid-Axford profiles is the 
same, and the only difference between them is the time dependence outside the 
exponential function. While the Reid-Axford profile has time dependence, the 
inverse Gaussian profile has time dependence. This difference is ascribed 

to the different physical mechanisms considered for the two profiles. While the 
Reid-Axford profile is derived on the basis of lateral diffusion in the solar corona, 
the inverse Gaussian is based on vertical drift accompanied by diffusion in the 
solar corona. In reality, accelerated particles can drift and diffuse anisotropically 
in three dimensions, and the degrees of drifting and diffusion along each direction 
depend on the coronal magnetic field configuration, therefore, the injection model 
that is more suitable depends on the event. 

As mentioned in the previous section, an identical SEP injection profile yields 
different intensity profiles depending on the scattering mean free path in inter¬ 
planetary space. The mean free path is a crucial parameter for determining the 
SEP transport conditions. However, the comparison shown in the previous section 
predicted almost identical SEP intensity prohles over a wide range of mean free 
paths. This is because, in the comparison the injection profile was adjusted to fit 
the observational data for each mean free path. In contrast, the anisotropy profiles 
are widely scattered for similar calculated intensity profiles and different mean 
free paths. Figure |9] illustrates the calculated anisotropy profiles for the event on 
27th January 2012. The lines from bottom to top represent the anisotropy profiles 
for radial mean free paths ranging from 0.2 to 2.0 AU with increments of 0.3 AU. 
The anisotropy varies considerably with the mean free path, although the calcu¬ 
lated intensity prohles are similar. Although no observation of the anisotropy has 
yet been reported for SEPs with energy greater than several tens of MeV, if such 
observations become available, our model may be able to estimate the scattering 
mean free path in interplanetary space in near real-time, which will be valuable 
for the study of SEP transport and for further predicting the SEP intensity. Thus, 
anisotropy observations for SEPs, particularly with energy greater than several 
tens of MeV, will be valuable for SEP research and prediction. 

Figure [To] illustrates the derived injection prohles for the event on 27th Jan¬ 
uary 2012. Golored lines represent the injection prohles for radial mean free paths 
ranging from 0.2 to 2.0 AU with increments of 0.3 AU. We can recognize from 
the hgure that the longer the mean free path, the longer the injection time scale. 
A long mean free path such as 2.0 AU results in a quite long injection time. One 
possible interpretation of the long injection time is that the SEP source is far from 
the magnetic foot point of the Earth. In this case, the SEPs must diffuse across 
the magnetic held lines in the solar corona to arrive at the magnetic held lines 
connected to the Earth. It takes the SEPs a long time to escape from the upper 
corona to interplanetary space. 

Because the simulation is spatially one-dimensional, perpendicular dihusion 
and gradient-curvature drift are not included, and particles only move along mag¬ 
netic held lines. This means that particle acceleration sites must be connected to 
observation sites via magnetic held lines, which may limit the applicability of the 



12 


Yuki Kubo et al. 


method to well-connected events. All events analyzed in the previous section were 
accompanied by solar flares in the western hemisphere of the Sun, suggesting that 
these events were well-connected. On the other hand, according to Gopalswamy et 
al. [2013], the event on 27th January 2012 was not well-connected because of the 
large latitudinal separation (~ 40°) between the CME location and the ecliptic 
plane, although this event occurred in the western hemisphere of the Sun. This is 
why the initial phase of the event is slightly slow-rising. As already mentioned, the 
event on 6th January 2014 occurred at the far-western side of the Sun, implying 
that this event was also not well-connected. Therefore, our simulation sometimes 
works well by choosing a slow-rising injection time profile even though the SEP 
events are not well-connected, as long as the accompanying flare/CME occurs in 
the western hemisphere of the Sun. 

Our model may not reproduce events with a plateau around the intensity peak, 
such as the event on 29th September 1989, which occurred behind the western limb. 
This is because these events often include interplanetary-shock-accelerated SEPs, 
while our model only deals with SEPs accelerated near the Sun (at a solar flare 
and/or in the solar corona). As an interplanetary shock will seldom accelerate SEPs 
up to several tens of MeV (although an extremely strong interplanetary shock can 
achieve this), most high-energy SEPs, such as those observed by the GOES P6 and 
P7 channels, are accelerated near the Sun, and the intensity time profiles observed 
by the GOES P6 and P7 channels do not show a plateau. Actually, while the GOES 
PI to P5 channel data for the 29th September 1989 event showed plateau prohles, 
the P6 and P7 channels showed decay profiles rather than plateau profiles. As this 
study focused on high-energy SEPs, such as those observed by the GOES P6 and 
P7 channels, our model will work well even for the event on 29th September 1989. 
Erom 1997 to 2014, 29 events in the western hemisphere (longitude greater than 
W30) with intensities of SEPs with energy greater than 100 MeV reaching 1 PFU 
were observed by the GOES satellite. By visually checking the GOES P6 and P7 
channel data for the 29 events, we concluded that 23 of the 29 events can be fitted 
using our model (almost 80% of the western hemisphere events). 

As an application of the numerical simulation, we may apply the calculated 
results to predict the decay phase time prohle of an SEP event. We fitted the results 
of the calculation to the observed initial phase data of the GOES P6 channel with 
radial mean free paths ranging from 0.2 to 2.0 AU with increments of 0.3 AU. The 
mean m and standard deviation a of the inverse Gaussian injection profile were 
adjusted to match the simulated intensity of the observed initial phase of the SEP 
event on 27th January 2012. After the fitting procedure, the decay phase intensity 
of the SEP event was calculated using the adjusted parameters m and a. The 
results are shown in Figure [11] The red and green points indicate the observation 
data used and those not used for fitting, respectively. Because the mean free path 
influences the predicted profile, the plotted predicted prohle (blue line) is the 
average of the prohles predicted at different mean free paths. The upper and lower 
envelopes of all predicted prohles can be regarded as the upper and lower limits 
of the prediction, respectively. Thus, the region between these limits (gray region) 
can be regarded as the error bars or error range of the prediction. Because most 
of the observed data not used for htting (green points) are located inside the error 
range in Figure (TT] our calculations reproduce the decay phase of this event with 
a small error range. This is one of the applications of our numerical simulation. 
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suggesting the possibility of using the SEP transport simulation to predict SEP 
decay phase time profiles. 

As mentioned in the previous section, our simulation can closely reproduce 
GOES P7 channel data. The energy of the P7 channel is just the energy that pro¬ 
duces an extreme event such as a GLE. Therefore, our simulation is expected to be 
useful for reproducing and/or predicting GLEs and radiation doses for astronauts 
and aircrews. 


6 Summary 

We developed a numerical simulation code to calculate SEP intensity profiles as one 
component of a recently developed warning system for aviation exposure to SEPs. 
The code solves a spatially one-dimensional focused transport equation recast as 
a set of stochastic differential equations. An inverse Gaussian distribution was 
employed instead of the conventionally used profile for the injection time profile 
of SEPs near the Sun. The code provides the temporal profile of the differential 
intensity of an SEP event. To validate the code, the simulation results obtained 
using the code were fitted to observational data for four SEP events on 27th 
January 2012, 13th March 2012, 17th May 2012, and 6th January 2014 that were 
recorded by the GOES P6 and P7 channels, the square error of the fitting was 
minimized using the simplex method. The model predicted all the event profiles 
with satisfactory accuracy, indicating that the inverse Gaussian distribution can 
be used as the SEP injection profile instead of the well-known injection profile. 
By successfully reproducing SEP events, this code is expected to be applied to 
develop the WASAVIES warning system for aviation exposure to solar energetic 
particles. 
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Fig. 1 Non-orthogonal curvilinear coordinate system (s, 0, (p) for describing the focused trans¬ 
port equation. The solid and dotted blue lines are respectively a Parker magnetic field and its 
projection onto the X-Y plane, s is the distance of P measured from the origin O along the 
magnetic field line, 9 is the polar angle at P measured from the Z axis, (f) is the phase angle of 
the magnetic field line at O measured from the X axis, which is different from the phase angle 
</? at P, r is the radial distance from O to P, e^, e^, and are the basis vectors at position 
P, and (r, 0, (p) constitutes the spherical coordinate system. 
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trum at the Sun. Colored lines show the temporal evolutions of spectra at the Earth. The 
spectra are normalized using a value at 1,000 MV. 
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Fig. 5 Comparison between calculated SEP intensity and the GOES P6 (upper panel) and 
P7 (lower panel) channel data for the event on 27th January 2012. The gray points are the 
observation data, which are plotted in the figure by thinning them out to increase their visibil¬ 
ity. The colored lines show the SEP intensity profiles calculated using our method for various 
mean free paths. 
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Fig. 6 Same as Figure [S] but for the SEP event occurring on 13th March 2012. 
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Fig. 7 Same as Figure [S] but for the SEP event occurring on 17th May 2012. 
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Fig. 8 Same as Figure (5] but for the SEP event occurring on 6th January 2014. 
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Fig. 9 Calculated anisotropy profiles for the event occurring on 27th January 2012. From 
bottom to top, the lines are plotted for mean free paths ranging from 0.2 to 2.0 AU in 0.3 AU 
increments. 
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Fig. 10 Derived injection profiles for the event occurring on 27th January 2012. The lines are 
plotted for mean free paths ranging from 0.2 to 2.0 AU in 0.3 AU increments. 
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Fig. 11 Results of GOES P6 channel data fitting using only the initial phase data of the 27th 
January 2012 event. The red and green points are the observation data used and those not 
used for the fitting, respectively. The blue line shows the SEP intensity profile calculated using 
our simulation, which was averaged over the results for mean free paths ranging from 0.2 to 
2.0 AU with increments of 0.3 AU. The gray area shows the error range (see text for details). 






